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We use a post-Newtonian diagnostic tool to examine numerically generated quasiequilibrium 
initial data sets for non-spinning double neutron star and neutron star-black hole binary systems. 
The PN equations include the eflFects of tidal interactions, parametrized by the compactness of the 
neutron stars and by suitable values of "apsidal" constants, which measure the degree of distortion 
of stars subjected to tidal forces. We find that the post-Newtonian diagnostic agrees well with the 
double neutron star initial data, typically to better than half a percent except where tidal distortions 
. . - are becoming extreme. We show that the differences could be interpreted as representing small 

residual eccentricity in the initial orbits. In comparing the diagnostic with preliminary numerical 
data on neutron star-black hole binaries, we find less agreement. 
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Advances in numerical general relativity have made it possible to solve Einstein's equations for the astrophysically 
important problem of the inspiral and merger of binary systems of compact bodies (black holes or neutron stars) 
with unprecedented accuracy and robustness, both in the initial data obtained from Einstein's equations, and in the 



qh| subsequent time evolution of the system. Work is ongoing to stitch together results from post-Newtonian theory, 

^ ■ which accurately describes the early part of the inspiral, with numerical relativity, which describes the final few 

{^JO orbits, merger and ringdown of the final black hole, in an effort to develop a complete description of the evolution and 

I— I , gravitational- wave emission from such systems that can be used effectively in analysis of data from gravitational- wave 
interferometers 0, [3, BI3]- 

^ , As part of this program to link post-Newtonian theory with numerical relativity, we proposed and developed a 

' diagnostic tool, based on the post-Newtonian (PN) approximation, for analysing numerical relativity initial data sets 

00 " . . . mmn „ . 



for compact binaries [E, S 0] ■ Recall that the post-Newtonian approximation is effectively an expansion of Einstein's 



■ theory in powers of a small parameter e ~ (w/c)^ ~ Gm/rc^. 
' Numerical data sets, obtained from solutions of the initial-value equations of general relativity, are generally config- 
. ured to represent the quasicircular orbit of a compact binary system late in its inspiral phase, when the bodies have 

■ only a few orbits remaining before merger. Such a quasicircular orbit (circular, apart from the radiation-induced in- 
, spiral) is the expected end-point of evolution of a compact binary system under the circularizing effect of gravitational 

■ radiation reaction in the absence of external perturbations. 

] The PN diagnostic provides analytic expressions for the total binding energy (the difference between the total 
' ^ ' gravitational mass of the system and that of the two stars in isolation) and angular momentum of the system as 
Kn] ■ a function of its orbital angular velocity, and allows for arbitrary mass ratios, spins, finite-size effects such as tidal 
^ interactions, and a non-zero orbital eccentricity. We found that the PN expressions for circular orbits agreed with the 
■ - - ' numerical relativity results surprisingly well (see also Q), even for orbital separations in a highly relativistic regime 
where e is not so small, but that there were systematic differences that could not be blithely attributed to errors in 
either the PN approximation or the numerical data. 

For binary black hole initial data, for example, we found [^,0] that our diagnostic gave a better fit to the numerical 
data from several groups if the orbital eccentricity were allowed to be non-zero, with values as large as 0.03. Here 
the bodies were viewed as residing initially at the apocenter of the orbit, in keeping with the constraint built into 
the numerical initial data that the bodies be moving perpendicularly to their separation. And indeed, subsequent 
calculations of the time evolution of these initial data using numerical relativity showed that the orbits were slightly 
eccentric, with values of eccentricity agreeing in order of magnitude with those predicted by the PN diagnostic (see, 
eg. 0). Efforts are now being made to "tweak" the initial data for binary black holes in order to achieve orbits that 
more nearly approximate the expected quasicircular orbits [lol . [ill , [l^ . 
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In this paper we turn our attention to binary systems containing non-rotating neutron stars (NS). In Idl we applied 
the diagnostic to one neutron star binary initial data set for corotating stars obtained by Miller et al. Il3|. but here 
we focus on an extensive set of non-spinning binary neutron star initial data obtained by Taniguchi et al. [l4l. [isL [l^ , 
including both equal-mass and unequal-mass systems; and on preliminary initial data sets for neutron star-black hole 
binaries obtained by Taniguchi et al. [lit . 1I81 . il9l] and Grandclement (20| . 

While for black holes, tidal effects are negligible at the separations in question, for neutron stars, tidally induced 
distortions must be taken into account (see [6| for discussion). These effects depend on the size of the neutron star, 
as encoded in a "compactness" factor, M/R, where M and R represent mass and radius, and on a set of "apsidal 
constants", ki, which measure the degree of distortion of the body as a result of an external tidal force, for each 
multipole order I. These are also known as "Love numbers" in other contexts. The apsidal constants depend on the 
equation of state (EOS) used to model the neutron star interiors. 

For given values of the masses, compactness parameters, apsidal constants and eccentricity, the PN expressions for 
binding energy and angular momentum can be evaluated and compared with the data. Alternatively, for example, 
one can leave the eccentricity as a free parameter and solve for that value that gives the best agreement between the 
PN and the numerical results for the given angular velocity. 





FIG. 1: Best-fit eccentricity estimated from the energy (thick lines) and the angular momentum (thin lines) for irrotational, 
equal-mass neutron star-neutron star binaries. Top: we fix F = 2.00 and consider different values of the compactness. Bottom: 
we fix M/R = 0.14 and consider different values of F. Results on the left use the Newtonian apsidal constants, results on the 
right use the "relativistically corrected" apsidal constants of Eq. (|A6|) . 



In comparing with the non-spinning NS-NS initial data sets of [3, [H, [iBl , our main results can be summarized 
using Fig. [TJ These figures show the eccentricity required to match PN theory with the numerical data for equal-mass 
systems, for both binding energy (thick lines) and angular momentum (thin lines), for various equations of state, as 
represented by the adiabatic index F, and for various values of the neutron star compactness M/R. The panels on 
the left use values for the apsidal constants k2 and ^3 calculated using Newtonian theory, while those on the right 
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use values for the dominant k2 apsidal constant estimated from rotating relativistic stellar models (see Appendix [X)) . 
The figures are plotted against mfi, where m is the total gravitational mass of the two isolated stars and is the 
orbital angular velocity. The left-hand end of each plot corresponds to the "Newtonian" limit, where the stars are 
farther apart; on the right-hand end, the sequences terminate when the stars are within 25 to 40 percent of touching, 
where cusps in the surface shape are beginning to develop, signalling the onset of mass shedding. 

As discussed below (Sec. HIl, the inferred eccentricity is also an approximate measure of the fractional difference 
between the data and PN theory assuming circular orbits, so the figures can equally well be interpreted as a measure 
of the accuracy of the PN-numerical comparison. 

From these figures, together with considerations in following sections, we summarize our main results. 

PN theory gives an excellent fit, better than half a percent in most cases, to the binding energy of the numerically 
generated initial data sequences for all but the closest separations. Tidal effects matter, and one must use the correct 
compactness factor and apsidal constants corresponding to the isolated stars and equations of state used in the 
numerical solutions. 

// circular orbits are assumed, there are systematic differences between the numerical data and the PN diagnostic. 
These could be interpreted as being the result of residual eccentricity in the simulations, albeit at levels smaller than 
those inferred from the BH-BH models in 0, 0|- The differences are larger than can be accounted for either by 
truncation errors in the PN series or by numerical errors in the initial data sets. 

When best-fit eccentricities are calculated, the values inferred from the energy diagnostic are generally smaller than 
those inferred from the angular momentum diagnostic. The same phenomenon was seen in the BH-BH diagnostics 
0] ■ One possibility is that there are poorly understood systematic effects in calculating total angular momentum from 
numerical initial data sets, or in fixing the intrinsic angular momenta of the individual stars to the desired values 
(zero in this case). 

A comparison with preliminary initial data for neutron star-black hole binaries shows that the required eccentricities 
are larger than for NS-NS binaries. Furthermore, they increase with the mass ratio tobh/wns- However, because 
these are early days for NS-BH simulations, the numerical errors in the models are larger than in the more mature 
NS-NS simulations, so these conclusions should be regarded as tentative. Our PN diagnostic could be a useful tool in 
helping to refine and improve these initial data sets. 

Because tidal effects fall off rapidly with distance and scale with the size of the bodies, the Newtonian tidal contribu- 
tions are already small, and thus we have argued that post-Newtonian tidal effects are not needed. However, because 
the neutron stars are highly relativistic, one might question the use of apsidal constants derived from Newtonian 
theory. We have addressed this by constructing models of isolated, slowly rotating, fully relativistic neutron stars, 
and by using the rotationally induced deformation to estimate the "relativistic" apsidal constant k2. The resulting 
values can vary by a factor of two over the range of compactness factors considered, leading to small but noticeable 
effects in the inferred eccentricities (right-hand panels of Fig. [T|) . 

The remainder of this paper provides details. In Sec. |TT] we summarize the post-Newtonian diagnostic equations 
derived in Q, apphcable to both NS-NS and NS-BH non-spinning binaries. In Sec. [TTIlwe consider the catalogue of 
quasiequilibrium NS-NS sequences c omp uted by the Meudon group [3, [l^, [l^ . Section IIVI considers NS-BH binaries 
obtained by Taniguchi et al. [TtI. [l^ IT9l| and Grandclement [2(|. Section |V] presents concluding remarks. Appendix El 
provides details of the calculation of relativistic corrections to the energy, angular momentum and quadrupole moment 
of rotating isolated neutron stars described by a polytropic equation of state. These calculations were performed using 
the Hartle-Thorne slow-rotation code described in [21], and are used to estimate "relativistically corrected" values of 
the apsidal constant /j2- Hereafter, we use units in which G = c = 1. 



We consider a binary system of bodies of mass toi and TO2, where these masses denote the total, or gravitational 
mass of an isolated, non-rotating body of given baryon number and equation of state. For black holes, the masses 
denote irreducible mass. We also define the total mass, m = mi -\- m2, and the reduced mass, ^ = mim2/m, along 
with 77 = fi/m — mim2/m^ (0 < 77 < 1/4). We define two quantities e and related to the eccentricity and semi-latus 
rectum of the orbit, according to: 



II. EQUATIONS OF THE PN DIAGNOSTIC 




(2.1) 
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where fip is the value of the orbital frequency Q where it passes through a local maximum (pericenter) , and fla is the 
value of f2 where it passes through the next local minimum (apocenter). These definitions are discussed and justified 
in [1,01; in the Newtonian limit, they agree exactly with the standard definitions. It follows from Eqs. (|2.ip that C 
also has the property that 



c 



2/3 



2/3 



(2.2) 



The diagnostic was restricted to systems of possibly rotating bodies with spin axes aligned perpendicular to the 
orbital plane (see 0, Eqs. (2.35), (3.6), (3.7), (3.15) and (3.17) for the complete set of ingredients). However, 
because we will be considering numerical initial data sets for systems with nonrotating bodies, the formulae simplify 
significantly, because there will be no contributions from rotational kinetic energy or angular momentum, spin-orbit 
or spin-spin coupling, effects due to rotational distortions or tidal-rotational coupling, and so on. As a result, the 
total binding energy E and total angular momentum ,/ of the system can be written 



E — -EHarm + £^Tidal , 
J = "/Harm + "^Tidal ■ 



(2.3) 



where the "point-mass" contributions to E and J, accurate to third post-Newtonian (3PN) order, and expressed in 
harmonic coordinates, are given by (recall that ( ^ 0(e)) 



-mri{l — e^)(' < 1 
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(2.4a) 



(2.4b) 



Expressions for E and J are also available in Arnowitt-Deser-Misner (ADM) coordinates "S] ; however the differences 
are negligible for the systems in question. The tidal contributions to E and J are given by. 



^^Tidal 
.^Tidal 



777.77(1 — e^) 
2 



18 



(9 + lOe^ - 3e^)BC^ + (13 + 496^ + 7e* - 5e^)CC^ 



24 



777,^77 



(3 + lOe^ + 3e4)i?C'/' + ^(1 + 76^ + 7e^ + e'')CC'^' 



(2.5a) 
(2.5b) 
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where 



B 
C 



6ri 
8t] 



V m / 



(2.6) 



For each body, Qa = Ra/nia is the inverse of the "compactness parameter", where Ra is its radius in harmonic 



coordinates, and fcj"'' ^'^d ^3"^ denote the "apsidal constants" for angular harmonics I — 2 and I — 3, respectively 
(higher harmonics make negligible contributions). In Newtonian gravity, the apsidal constants depend only on the 
harmonic index I and on the equation of state; their values are shown in Appendix |^ Table IIIII Notice that the 
leading contribution to the tidal terms is of order relative to the Newtonian contributions to E or J, or effectively 
5PN order. This is why a Newtonian treatment of tides is justified. We also note that, in calculating E/m and J/rri^, 
only mass ratios are relevant, since, apart from the dependence on C, the expressions depend only on 7/, mi/m or 
1712 /m. 

In comparing this analytical diagnostic with numerical initial data, one must carefully translate the meaning of the 
variables in the two approaches. The PN masses rua are the physically measured masses of the corresponding isolated 
body; these correspond directly to the ADM masses of the isolated bodies defined in numerical relativity. Similarly, 
numerical stellar radii are generally computed in ADM coordinates, which, for an isolated body corresponds to the 
isotropic radial coordinate of the Schwarzschild geometry; we denote this radius by Radm- Alternatively, one could 
compute, from the numerical data, the proper circumferential or areal radius of the isolated star [C/2i: or (^/47r)^/^], 
which would correspond to the Schwarzschild radial coordinate Rs- The relationship between the PN harmonic radial 
coordinate Rh and the others is 



R 



H 



Rs - M, 



ADM 



Ra 



DM 



^^ADM 



(2.7) 



Since for neutron stars, Madm / Radm < 0.2, the difference between Rh and Radm is generally less than one percent, 
and is thus negligible for our purposes. For non-rotating black holes, the relevant radius is Rh ~ "^bh, but since tidal 
effects are negligible for the separations to be considered, this will not play a role. 

Finally, we note that, inserting Eq. (|2.2p into the Newtonian part of Eqs. (|2.4p . we obtain to leading order. 



E/m w -^77(771^2)^/3 
J/7772 _ riim^)-^'^ 



l + _e + 0(e2) 
l-le + 0{e') 



(2.8) 



Thus, apart from factors of 4/3 or 2/3, the eccentricity e can be viewed as a measure of the fractional deviation of 
the energy and angular momentum from its circular orbit value for a given 777r2. If the inferred value of e is positive, 
the binary sytems is at apocenter; if the inferred value is negative, it really corresponds to a positive eccentricity, but 
with the binary system at pericenter. 



III. QUASIEQUILIBRIUM BINARY NEUTRON STAR INITIAL CONFIGURATIONS 

A. Numerical Calculations 



The first numerical calculations of initial data for binary neutron stars in full general relativity were carried out, 
for the case of momentarily corotating bodies, by Baumgarte et al. [^s', '2^1 and Marronetti et al. [23|. Miller et al. 
[is] performed convergence tests and included error bars in their simulations of corotating neutron stars. 

The first results for the more physically realistic case of non-spinning, or irrotational binary neutron stars were 
presented by the Meudon group using spectral methods [3, HI] , by Marronetti et al. using a single-domain finite 
difference method in Cartesian coordinates [29,], and by Uryu et al. using quasispherical coordinates [30, 31, 32]. 
Generalization to other values of the stellar rotation rates was carried out in [33] . 

We will focus on simulations of non-spinning binary neutron stars by the Meudon group in [ij, [l^ [l^. They 
assume that the system is in a quasiequilibrium state, meaning that, in a reference frame that momentarily rotates 
with the orbital angular velocity, the system is stationary (technically implying the existence of a "helicoidal Killing 
vector" d/dt + rtd/dcf)). They assume that the fluid flow is irrotational with respect to the global frame and that the 
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induced spatial metric is conformally flat. They use a multidomain spectral method with two patches of surface- fitting 
coordinates, to solve for the five metric functions (lapse, conformal factor and shift 3- vector) and the velocity potential 
of the fluid. The input information is an equation of state of each star; an initial coordinate separation d; and the 
central enthalpy (or density) of each star, or alternatively the baryonic mass of each star. 

The equation of state used is a zero-temperature polytrope with p — np^ , where where p and p are the pressure 
and mass density, and T is the adiabatic index. Such equations of state are also parametrized by a "polytropic index" 
n, related to F by F = 1 + 1/n. Most groups use "polytropic" units, whereby one defines a length scale 

= ^V[2(r-i)] ^ (3,1) 

and rescales every quantity by i?poiy to make it dimcnsionless. With this choice, equilibrium models are characterized 
only by F and by the (dimcnsionless) central energy density of each body (see eg. Sec. 2.6.1 of 34] for a detailed 
explanation of these units). 

In Ref. 14] calculations were limited to binaries of equal mass and to a single polytropic EOS with F = 2, 
M/R = 0.14 and baryonic mass Mb = 1.625Mq. Results were given in their Table IIF Taniguchi and Gourgoulhon 
[la ], again with F = 2, extended the calculations to corotating and non-spinning binaries of varying compactness, and 
also with unequal masses. Finally in Ref. [l^ calculations were presented for an extensive sample of values of F and 
of M/R, for both equal and unequal masses. Comparison of the endpoints of their sequences of models with those 
of Uryu et al. [13, HH, [s^ showed good agreement [lB| . Our Tables [T[ and [TT| summarize schematically the different 
quasiequilibrium sequences considered by the Meudon group, pointing to the relevant Tables in the original papers 

TABLE I: Filled entries mark combinations of F = and M/R considered by the Meudon group for equal-mass binaries. 

For non-spinning configurations, bullets mark cases for which we plotted E/m and J/rr? , and circles mark cases for which 
we plotted the inferred eccentricities. The third and fourth columns give the maximum baryonic mass M'^'^^ and maximum 
compactness (M/i?)™^'', respectively, for the given polytropic model. "Cor" means corotating stars, "NoS" means non-spinning 
stars. The last three columns give the Tables and References where the corresponding numerical data are tabulated. 



Parameters 


M/R 


Table 


Ref. 


r n MT" (M/i?)"""" 


0.08 0.10 0.12 0.14 0.16 0.18 0.20 


Cor NoS 




1.80 5/4 0.217 0.172 
2.00 1 0.180 0.214 
2.25 4/5 0.162 0.252 
2.50 2/3 0.155 0.278 


• X X •o 

• O #0 

X O X X 

• O X X • 


X XII 

I III 
VI VIII 

II IV 


fl6] 
fl5] 
fl6] 
fl6] 



TABLE II: Same as Table |TJ but for unequal-mass binaries. 



Parameters 


{M/R)i + {M/R)2 


Table 


Ref. 


r n 


(0.08 + 0.10) (0.10 + 0.12) (0.12 + 0.14) (0.14 + 0.16) (0.16 + 0.18) (0.18 + 0.20) 


Cor NoS 




1.80 5/4 
2.00 1 
2.25 4/5 
2.50 2/3 


XXX 

• X • 

XXX 

XXX 


XI XIII 

II IV 
VII IX 

III V 


fl6] 
[15] 
fl6] 
[16] 



B. PN diagnosis of Meudon data 

The tables of Refs. [H, [IBj display the total ADM mass M, angular momentum J and angular velocity f2 
of the system (among other variables) for the various models. Overbars denote that all variables are in polytropic 
units. One piece of information missing from those tables is the ADM mass of each star in isolation; these values were 
kindly provided to us by Keisuke Taniguchi, along with tables of all data that went into those references, but with 
full untruncated precision. 

It is then straightforward to compute E/m and J/rri^ using Eqs. (|2.3p - (|2.6p as functions of mfi, for each value of 
F (which determines the apsidal constants) and the compactness parameter, together with the provisional assumption 
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that e = 0. From the numerical data the corresponding binding energy and angular momentum are given by 



E 
m 
J 



M 



Num 



Num 
(mri)Num 



A^ADM,0 
J 

-"'-'adm,o 



(3.2) 



where Madm,o is the sum of the ADM masses of the two stars in isolation. We then plot (i?/TO)Num and ( J/m^)Num 
against (TOil)Num- Since the final quantities E/m, J/m^ and mil, are dimensionless, the polytropic units scale out. 
The compactness parameter is quoted as (AI/R), where here M is the ADM mass of the isolated star, and R is the 
areal, or Schwarzschild radius (the value is independent of whether polytropic or normal units are used); as a result, 
using Eq. (|2.7p , we can read off the value 



(9) 
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FIG. 2: Binding energy and angular momentum vs. mQ. for equal masses, F = 2.0, {M/R)i = {M/R)2 = 0.12. Top: circles 
are the numerical data points, solid (black) line is the PN diagnostic using the Newtonian apsidal constant values (Table 
IIII|) for r = 2, dashed (red) line has no tidal contributions; dot-dashed (blue) line uses the uniform density values for the 
apsidal constants; dotted (black) line uses the "relativistically corrected" apsidal constants of Eq. (|A6|l . Bottom: Numerical 
data (circles), point-mass (dashed-red) and uniform-density (dot-dashed-blue) plotted relative to F = 2 PN diagnostic as the 
baseline. Triangles denote numerical data plotted relative to F = 2 PN diagnostic using the relativistically corrected apsidal 
constants. 
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FIG. 3: Same as Fig. H for F = 2.0, (M/ii)i = {M/R)2 = 0.18. 



Figures [2] and [3] show the results for equal masses, F = 2 and compactness values 0.12 and 0.18, corresponding 
respectively to neutron star harmonic radii 7.3 and 4.5 in units of the stellar mass. The top figures show, on a gross 
scale, the overall agreement between the PN and numerical results. Using differences between data sets, the bottom 
figures magnify the scale: the circles and triangles show the difference between the numerical data and our diagnostic 
using the Newtonian and relativistically corrected apsidal constants respectively. The dashed and dot-dashed curves 
illustrate the importance of taking tidal effects into account; respectively, they show the difference between the PN 
values for energy and angular momentum assuming the Newtonian apsidal constants and those assuming point masses 
(no tidal effects) or those assuming uniform density bodies (maximum tidal effects) . It is clear that the numerical data 
agree much better with PN formulae that include appropriate tidal effects than they do with formulae that assume 
either point masses or homogeneous bodies. The agreement is at the level of a percent or better over the range of 
mJ7 shown. However, as in the binary BH case, there appears to be a systematic offset between the data and the 
diagnostic for J /m? . We note that the use of a relativistically corrected apsidal constant (triangles) has a small effect 
on the agreement. As discussed in Appendix [Xl relativistic effects tend to decrease the apsidal constant k2', this has 
the effect of raising the data points on the difference plots compared to the Newtonian diagnostic baseline. 

Figure[4]shows differences relative to the PN diagnostic baseline for other values of F and M / R for equal masses. The 
agreement between diagnostic and data seems to worsen as the EOS progress from soft to hard and the compactness 
increases; again there is a systematic offset in J/m^, even in the non-relativistic limit. Figure [H] shows differences for 
selected configurations of unequal-mass neutron stars. 

Instead of comparing the data with a circular-orbit diagnostic, we now let the eccentricity be a free parameter, and 
find the value of e that gives the best fit to the data. Using both the Newtonian and the relativistic apsidal constants, 
we show the results for a range of EOS and compactness parameters in Fig. [TJ We note that the inferred eccentricities 
are generally smaller than those inferred from the binary BH initial data, and that there is again a difference between 
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FIG. 4; Same as the bottom panels of Fig [21 for various values of F and M/ R, for equal mass binaries. 
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FIG. 5: Same as Fig. |4j but for binaries with different masses. 



the values inferred from angular momentum and those inferred from energy. The energy-inferred eccentricities are 
generally less than 0.005 for all but the high-frequency ends of each quasiequilibrium sequence. The end points of 
these sequences correspond to binary separations only 25 to 40 percent larger than the sum of the stellar radii; at 
these separations, tidal effects are sufficiently strong that cusps in the surface shape of the stars are about to form, 
signalling the onset of mass shedding. Not only is our simple tidal model likely to be inadequate, but also there are 
larger errors in the numerical data near these end points (see for discussion). 

The inferred eccentricities are very small, and it is worth asking how significant they are, given that the PN 
approximation used in our diagnostic has truncation errors, and that the numerical simulations have numerical errors. 
Figure [S] is an attempt to quantify this. We plot the energy-inferred eccentricities for two representative models, 
with r = 2, {M/R) 0.12, and T = 2.5, {M/R) = 0.14. The solid curve is ten times what one would crudely 
estimate for the 4PN contributions to the energy for equal masses, or 10 X CV8 = 1.25(7^17)10/3^ while the dotted 
curves are (3/64)(7^fc2(mri)i'*/'^ for the two corresponding values of q and A:2, which represent a post-Newtonian- 
tidal correction, or Q times the leading tidal term in Eq. (|2.5ap . Also plotted as dashed curves with data points 
are the corresponding "virial" errors VE{M), computed from the numerical data, and tabulated in [isl. [l6j. where 
VE{M) = |Madm — -^Komarl/MADMl, whcrc MKomar IS a Certain global integral known as the Komar mass. This 
and similar virial theorems are one indication of the accuracy of the numerical solution. We see that PN truncation 
errors are smaller by an order of magnitude than the inferred eccentricities or "errors" from the diagnostic, and that 
the numerical errors are generally smaller, except in the extreme Newtonian limit. This gives some confidence that 
our PN diagnostic can in many cases "diagnose" real physics in the numerical simulations. 
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FIG. 6: Comparison of PN diagnostic eccentricities or errors with truncation errors in the PN approximation and with "virial" 
errors in the numerical simulations. 



IV. QUASIEQUILIBRIUM CALCULATIONS OF NEUTRON STAR-BLACK HOLE SYSTEMS 



The first numerical calculations of fully relativistic neutron star-black hole binaries were carried out by the Illinois 
group, for corotating stars, and in the limit where the black hole was 10 times more massive than the neutron star ^35j . 
However, following a number of technical improvements, the group recently presented results for quasiequilibrium 
sequences of non-spinning NS-BH binaries with a range of five mass ratios (including equal masses) , and two values 
of the compactness parameter of the neutron star jl7l . [l8| . They assume a conformally flat spatial metric and 
solve the initial value equations for the metric variables using multidomain spectral methods. Since then, they have 
also implemented a more accurate condition for ensuring that the black hole is nonspinning, first introduced in 
Ref. [36|. Tani guc hi has provided us with preliminary data for the improved sequences [19|. Using similar methods, 
Grandclement 20] also obtained quasiequilibrium NS-BH sequences, for a mass ratio of five and for four values of the 
compactness parameter. To date, all NS-BH calculations have assumed F = 2. 

As in the NS-NS case, results are given in polytropic units, with E/m, J/rri^ and mH. already tabulated as in Eq. (|3.2p 
along with the mass ratio M^^/AI^^^^ q and the ADM mass M^pjyj q and compactness parameter M^^j^ o/Ro of the 
isolated neutron star, where Rq is again the areal radius. From these data, we compute the factor q for the neutron 
star from q = i?oA^ADM,o ~ 1- 
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FIG. 7: Binding energy and angular momentum vs. mf2 for a black hole-neutron star binary with mass ratio 5. Plotted are 
data from Taniguchi et al. (circles) and from Grandclement (triangles), along with the corresponding PN diagnostic curve for 
circular orbits. 



Figure[7]shows the gross agreement between the Taniguchi et al. and Grandclement data, for the one case (mass ratio 
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5, compactness ~ 0.15) where they can be compared, and between the data and the corresponding PN diagnostics. 
Figure [5] shows the inferred eccentricities (or the accuracies of the PN comparisons) for all the NS-BH data sets. The 
eccentricities are somewhat larger than in the NS-NS case (and significantly larger for the Grandclement mass ratio 
5 : 1 data), and vary more strongly as a function of mfi, especially near the end-points of the sequences, where the 
larger tidal deformations cause numerical errors. In addition, the eccentricities (or errors) appear to increase with 
mass ratio. This may not be surprising, as both numerical groups acknowledge that numerical errors tend to increase 
with mass ratio. At the same time, the convergence of the PN series depends on mass ratio, being best for equal 
masses, and rather poor for extreme mass ratios (see Fig. 1 of Ref. ^6]). In the equal-mass case, where the numerical 
errors are smallest and the PN series converges the best, the inferred eccentricities are small, less than 0.01 from the 
energy. 
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FIG. 8: Best-fit eccentricity estimated from the energy (thick lines) and the angular momentum (thin lines) for irrotational, 
black hole-neutron star binaries. Left: improved data from Taniguchi et al. 19j, for a compactness factor of 0.145, and four 
mass ratios. Right: data from Grandclement, for a single mass ratio of 5:1, and a range of compactness factors. The dotted 
curves in each figure correspond to a mass ratio 5 and similar compactness factors; they also correspond to the data plotted in 

Fig.m 



Because these NS-BH numerical simulations are still in their infancy compared to the BH-BH and NS-NS simula- 
tions, it is premature to draw definite conclusions from comparisons with PN theory. Instead, a PN diagnostic may 
be a useful tool in testing improved simulations. For example, in cornparing eccentricities inferred from the original 
data published in with those inferred from the unpublished data [loj that used an improved black-hole zero spin 
condition (plotted in the left panel of Fig. [5]) we found a significant reduction in the disagreement between values of e 
inferred from E and those inferred from J, especially at the closer separations. In Fig. [5] (left panel) the eccentricities 
from E and J closely track each other for the mass ratios 1, 2 and 3, both increasing slightly with mil; with the 
earlier data, the eccentricities inferred from J actually grew strongly with mfi, diverging from those inferred from 
E. A similar improvement resulting from the new black- hole spin condition was noted in BH-BH simulations (see 0], 
Fig. 2). 



V. CONCLUSIONS 



We have used a post-Newtonian diagnostic tool to examine initial data sets for non-spinning, double neutron star 
and neutron star-black hole binary systems. The PN equations included the effects of tidal interactions, parametrized 
by the compactness of the neutron stars and by suitable values of apsidal constants. 

We found that PN theory agrees quite well with the NS-NS initial data, typically to better than half a percent except 
where tidal distortions are becoming extreme. The differences between PN theory and the numerics are sufficiently 
large and systematic that they could be interpreted as representing residual eccentricity in the initial orbits. It will 
be interesting to see if simulations of the time-evolution from these initial data sets show signs of eccentricity in the 
orbits. 

Comparison of our PN diagnostic with NS-BH initial data sets showed poorer agreement, not inconsistent with the 
larger numerical errors present in these preliminary results. We argued that a PN diagnostic could be a useful tool 
for examining future improved data for NS-BH binaries. 
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Tidal effects depend on the values of the apsidal constants, and we addressed how relativistic neutron-star structure 
would affect their values. An open question is whether our estimate of a relativistic k2 from rotational deformations 
is a reasonable estimate of the tidally- induced /c2 . This would be an interesting and important problem in relativistic 
stellar structure. 
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APPENDIX A: TIDAL EFFECTS IN RELATIVISTIC NEUTRON STARS 

In we argued that, for compact bodies such as neutron stars or black holes, the effects of finite size, such as 
rotational kinetic energy, tidal distortions, and so on, would be of effectively higher PN order than expected a priori, 
because the bodies' sizes scale as their masses M, i.e. R ^ q M, where g is a factor of order unity. Thus, for example, 
the rotational kinetic energy of a body can be expressed as i?Rot ^ Ilu'^/2 ^ MB?{M /r^){uj /Vt)"^ < E^q^iM /r)'^ , 
where Em is the Newtonian orbital energy, and w is the body's rotational angular velocity. For compact bodies that 
are rotating no faster than the orbital angular velocity, rotational kinetic energy can thus be viewed as an effectively 
2PN contribution to the energy, even though formally it is a Newtonian effect. 

Similarly, tidal deformations have the dependence i?Tidai ~ (M^/i?)(i?/r)^ ~ ENq^{M/r)'^ , and so are effectively 
5PN order. However, for neutron stars, q can be as large as 8, and so the effects could be important, as we have seen. 

While we have argued that the smallness of these effects justifies a Newtonian treatment of tidal and rotational 
deformations f6l|, one might worry that, in relativistic neutron stars, the values of such quantities as apsidal constants 
could differ strongly from their Newtonian values. In this appendix, we address this question. 

Unfortunately, it is difficult to treat tidal effects in relativistic stars, because that would involve constructing 
a relativistic stellar model in the gravitational field of another body, which is exactly the problem of inspiralling 
binaries treated by numerical relativity. 

On the other hand, there is now a rather complete description of isolated rotating relativistic stars. In Newtonian 
gravity, there is no distinction between rotational flattening and tidal distortion, at least for I = 2 linear perturbations, 
and thus a single apsidal constant fc2 applies to both distortions (the only difference is that one distortion is oblate 
while the other is prolate). For higher I this degeneracy breaks down. In relativity, by contrast, there is a difference 
in principle between rotational and tidal deformation for 1 — 2, even for linear perturbations, because of the presence 
in rotating stars of such phenomena as frame-dragging and the relativistic contributions of rotational kinetic energy 
to self-gravity. Nevertheless, for the slowly rotating configurations that we expect to be relevant to estimating 
apsidal constants for small distortions, we might assume that the k2 that governs rotationally induced distortions in 
a relativistic star might not be too different from the k2 that governs tidal distortions. This permits one to read off 
k2 from the quadrupole moment induced by the slow rotation. 

Because the effects of fca were extremely small in our diagnostic, we will not consider this apsidal constant further. 
Therefore, in this appendix we will endeavor to estimate k2 for relativistic neutron stars and compare the results with 
those derived from Newtonian gravity. First we review the calculation of Newtonian apsidal constants. 

1. Apsidal constants in Newtonian gravity 

In Newtonian gravity, apsidal constants can be computed by a standard method for different polytropic indices n 
or r = 1 + 1/n and different values of I using the numerical method described in |22]. The results are listed, for 
example, in the classical monographs by Kopal [H, |2j| . Unfortunately these works do not list apsidal constants for 
n = 2/3, 0.8 and 1.25. Since these values are needed for neutron stars, we report here results of our own numerical 
calculations. In Newtonian gravity, apsidal constants are defined by the relation (see Appendix B of Q for a detailed 
discussion) 



' 2l + 2iji{R) 



(Al) 
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TABLE III: Newtonian apsidal constants (compare with Table II in j22j| and Table II in j23|): I is the angular harmonic index, 
r and n are adiabatic indices, with F = 1 + 1 /n. 



1 n = 2/3 n = 0.80 n = 1.00 n = 1.25 n = 1.50 n = 2.00 n = 2.50 

r = 2.5 r = 2.25 r = 2 r = i.8 r = i.67 r = i.5 r = i.4 

2 0.375966 0.325098 0.259909 0.194339 0.143279 0.073938 0.034852 

3 0.164696 0.138660 0.106454 0.075590 0.052849 0.024394 0.010192 

4 0.098546 0.081155 0.060241 0.040967 0.027393 0.011508 0.004342 

5 0.067424 0.054485 0.039293 0.025748 0.016569 0.006420 0.002220 

6 0.049791 0.039575 0.027827 0.017651 0.010984 0.003966 0.001272 

7 0.038649 0.030270 0.020810 0.012824 0.007745 0.002628 0.000789 



where R is the stellar radius and the function rji is a solution of the Clairaut-Radau differential equation 

rv'i + mim - 1) + QD{rji + 1) = l{l + 1) , (A2) 

with initial condition rii{0) — I ~ 2. A prime denotes a derivative with respect to the distance r from the center of the 
star; the quantity D = p{r)/p{r), where p{r) is the density of the undistorted configuration at r and p{r) is the mean 
density inside the volume of radius r. Of course, to obtain the density profile we must specify an equation of state. 
The simplest choice is a (Newtonian) polytropic EOS p = PcJ/", where pc is the central density, n is the polytropic 
index and y satisfies the Lane-Emden differential equation 

y" + -2/' + y" =0, (A3) 

r 

subject to the initial conditions j/(0) — I, y' (0) = 0. We integrated the system (jA2IIA3|) following the method described 
in [22|. Table Hill shows the results. When a comparison with Table II of [l^l is possible, the agreement is usually at 
the level of five to six decimal places or better. 

2. Apsidal constant k2 for relativistic stars 

We now turn to rotating relativistic stars, with the goal of usingtheir rotational distortions to estimate ^2. We use 
the general relativistic, slow-rotation Hartle-Thorne framework ^^iHIIj as described in f2l|. Rotating stars can be 
characterized by a dimensionless parameter e, given by 

e^LUy/W/M, (A4) 

which is the ratio between the body's angular velocity and the Keplerian angular velocity at the equator of a non- 
rotating stellar model (this e should not be confused with the post-Newtonian expansion parameter). Even for 
the fastest-rotating configurations (that is, when we consider corotating binaries at the ISCO) e ^ 0.35 (see the last 
column of Table IV)) . so a slow-rotation approximation is valid for our purposes. For constant-baryonic mass sequences, 
the Hartle-Thorne slow-rotation approximation induces an error with respect to the "true" value of the quadrupole 
moment which depends on the stellar mass and on the EOS. Since the apsidal constant ^2 is proportional to the 
star's quadrupole moment Q (see below)jOur calculation is valid as long as the error in Q induced by the slow- 
rotation approximation is small. In Ref. |2ll | we compared the Hartle-Thorne approximation with rotating neutron 
star calculations in full general relativity (see eg. Fig. 1 in 21]). From that comparison we conclude that our Hartle- 
Thorne based values of Q are good to within a few percent of the "true" values, which is more than acceptable for 
our estimates. 

Following standard conventions (see eg. (3^). we use the relativistic polytropic EOS p — np^ , with energy density 
e = p-\- P/{T — 1), and work in polytropic units. For n ~ 0.5 — 1 one obtains models with bulk properties comparable 
with those of observed neutron stars; polytropic indices n < 1 (n > 1) yield stiff (soft) stellar models, respectively. 
We are interested in fitting numerical sequences of quasiequilibrium binary systems containing neutron stars, which 
are computed at constant baryonic mass. Correspondingly, we consider rotational corrections along sequences of 
(isolated) rotating stars with constant baryonic mass. 

Given a sequence of rotating models, we can compute the "relativistic" apsidal constant fc2, defined as 
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for fixed values of M / R, where Q is the quadrupole moment of the star as determined from the exterior geometry. 
The last step makes use of the fact that, in the Hartle-Thorne formalism, Q and all rotation-dependent quantities 
(scalingwith some power of uj) are computed at the "reference" angular velocity lo = y^M/R^. As shown in Appendix 
B.3 of 6], this is the most natural relativistic generalization of the Newtonian apsidal constant. In fact, we verified 
that in the Newtonian limit M/R the Hartle-Thorne code reproduces the Newtonian apsidal constants listed in 
Table HIl Notice that since Q (x e^, and k2 oc Q/u'^, k2 does not depend on the rotation rate, but only on M/R. A 
concise way of presenting the data is by a polynomial least-squares fit (basically a post- Newtonian series) of the form 

2 

yt2 = ^fc(^)(Af/i?)^ (A6) 

j=o 

The coefficients fcj"''' are listed in Table IIVI We also report the percentage errors of the fits with respect to the 
numerical calculations, Ak™^^ — max[(fc2 — fc2"™)/fc2"™]: the fit is in excellent agreement with the numerical data 
for all values of n and in the whole range of M/R. For this reason, we used the analytical fits when computing the 
relativistically corrected apsidal constants to use in the diagnostic equations. 

Over the range of compactness factors shown, fc2 decreases by a factor of order two. This gives a small but 
noticeable shift in the eccentricities inferred from the numerical data, as displayed in Figs. [2] and [3l because the 
diagnostic baseline shifts toward a point-mass model (smaller /C2) the data referred to that baseline (the triangles) 
shift upward. 

TABLE IV: Fitting coefficients for the "relativistic" apsidal constants as defined in Eq. ()A6|I . Tlie fifth column shows the 
maximum percentage error of the fit with respect to the numerical data, and the last column shows the maximum value of 
M/R used in the fit. 



n kf^ 


..(1) 


1.(2) 


AfcS^"''(%) (M/i?)""" 


1.25 0.194339 


-0.773249 


0.733840 


0.32 


0.172 


1.00 0.259909 


-0.956592 


1.052813 


0.21 


0.214 


0.80 0.325098 


-1.144687 


1.377584 


0.23 


0.252 


2/3 0.375966 


-1.296016 


1.634947 


-0.80 


0.278 



3. Mass and moment of inertia of relativistic stars 

For future applications of the PN diagnostic to binaries with rotating neutron stars, it will be useful to have 
estimates of the variation in the mass and angular momentum of isolated rotating neutron stars, as a function of 
compactness and angular velocity. 

We used the code described in [2l| to compute rotational corrections to the gravitational mass of the isolated 
neutron star as a function of rotation for different values of n. We assume that the non-rotating model has one of the 
compactness M/R values chosen in TableHl this fixes the value of the baryonic mass for our nonrotating model. Then 
we increase the rotation rate at fixed baryonic mass (along an "evolutionary sequence" ) and compute the rotationally 
corrected gravitational mass 

Mrot = (l+/3)Mnonrot. (A7) 

More details on our numerical procedure are given in Sec. 3.1 of plj (the main difference being that there we fix the 
angular momentum and gravitational mass, here we fix angular velocity and baryonic mass). 

For each value of n and of the nonrotating star's compactness {M/R), we compute 11 models (the nonrotating 
model and 10 models equispaced in the angular velocity uj). Then we perform a polynomial least-squares fit of /3 as 
a function of the star's angular velocity u) (expressed in polytropic units): 

3 

/3 = ^/?('=)w'=. (A8) 

k=l 



The fitting coefficients 13^''^ are listed in Table IVl 
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TABLE V: Fitting coefficients for f3{uj) as defined in Eqs. HA7|) and (|A8|I . {M/R)c=o is the compactness and M is the 
gravitational mass (in polytropic units) of the nonrotating model. The seventh column shows the maximum percentage error 
of the fit with respect to the numerical data, and the last two columns show the maximum value of or e used in the fit. 



n 


(M/i?),=o 


M 


10=^ X 








) tD"^"" 


^max 


1.25 


0.08 


0.143165 


1.500742 


0.321866 


1.617628 


0.98 


0.048 


0.316 


1.25 


0.10 


0.164915 


0.138509 


0.325523 


0.740812 


0.49 


0.058 


0.316 


1.25 


0.12 


0.182153 


0.277598 


0.272061 


0.736909 


0.46 


0.071 


0.327 


1.25 


0.14 


0.194658 


1.062211 


0.207540 


0.858382 


0.78 


0.083 


0.329 


1.00 


0.12 


0.122286 


-0.072168 


0.153656 


0.082648 


2.61 


0.109 


0.331 


1.00 


0.14 


0.136233 


0.628680 


0.126834 


0.168151 


0.31 


0.125 


0.339 


1.00 


0.16 


0.147838 


0.781091 


0.113525 


0.172454 


1.69 


0.143 


0.347 


1.00 


0.18 


0.156736 


0.518670 


0.105391 


0.156755 


1.32 


0.160 


0.350 


0.80 


0.12 


0.087240 


-0.221682 


0.087986 


0.013138 


0.78 


0.152 


0.327 


0.80 


0.14 


0.099987 


0.534238 


0.076339 


0.051296 


-0.57 


0.171 


0.335 


0.80 


0.16 


0.111652 


0.278734 


0.078449 


0.034777 


-0.53 


0.191 


0.343 


0.80 


0.18 


0.122021 


0.345399 


0.075647 


0.035287 


-0.24 


0.208 


0.345 


2/3 


0.14 


0.080637 


0.794037 


0.049246 


0.034088 


0.59 


0.211 


0.332 


2/3 


0.16 


0.091591 


-0.248756 


0.061873 


-0.000499 


-0.37 


0.232 


0.340 


2/3 


0.18 


0.101844 


0.009904 


0.059124 


0.008106 


-0.43 


0.253 


0.347 


2/3 


0.20 


0.111209 


0.370593 


0.055858 


0.014234 


-0.15 


0.274 


0.351 



Relativistic corrections to the Newtonian moment of inertia, may be useful for the diagnosis of binaries containing 
rotating bodies, in order to evaluate the total angular momentum (including that of the stars themselves) and such 
effects as spin-orbit coupling. Following Q we define a dimensionless factor 

a = I/{MR^) = [Rg/Rf . (A9) 

The "radius of gyration" is defined as Rg ~ (//Af )^/^. In the Hartle-Thorne formalism, to leading order, Rg does not 
depend on the rotation rate (see the discussion following Eq. (11) of [^Tl). In Table IVTl we provide fitting coefficients 
for a post-Newtonian inspired power-series expansion for a, as computed using the Hartle-Thorne formalism: 

2 

a^^a^^\M/Ry . (AlO) 



TABLE VL Fitting coefficients for a{M/R) as defined in Eqs. (|A9|I and (lAlOp . The fifth column shows the maximum 
percentage error of the fit with respect to the numerical data, and the last column shows the maximum value of M / R used in 
the fit. 



Aq'"""(%) {M/Rf 



1.25 0.231893 0.256949 0.065809 0.19 0.172 

1.00 0.261297 0.274488 0.349978 0.15 0.214 

0.80 0.286378 0.275807 0.569004 -0.07 0.252 

2/3 0.303974 0.269094 0.714660 -0.19 0.278 
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